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We study the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction between two contact magnetic 
impurities placed on bilayer graphene (BLG). We compute the interaction mediated by the carriers 
of the pristine and biased BLG as well as the conduction electrons of the doped system. The 
results are obtained from the linear-response expression for the susceptibility written in terms of 
the integral over lattice Green's functions. For the unbiased system, we obtain some analytical 
expressions in terms of the Meijer G-functions, which consist of the product of two oscillatory terms, 
one coming from the interference between the two Dirac points and the second coming from the 
Fermi momentum. In particular, for the undoped BLG, the system exhibits the RKKY interaction 
commensurate with its bipartite nature as expected from the particle-hole symmetry of the system. 
Furthermore, we explore a beating pattern of oscillations of the RKKY interaction in a highly doped 
BLG system within the four-band continuum model. Besides, we discuss the discrepancy between 
the short-range RKKY interaction calculated from the two-band model and that obtained from the 
four-band continuum model. The final results for the applied gate voltage are obtained numerically 
and are fitted with the functional forms based on the results for the unbiased case. In this case, we 
show that the long-range behavior is scaled with a momentum that depends on Fermi energy and 
gate voltage, allowing the possibility of tuning of the RKKY interaction by gate voltage. 

PACS numbers: 81.05.ue; 75.30.Hx; 75.78.-n 



I. INTRODUCTION 

Graphene, a two-dimensional (2D) honeycomb lat- 
tice of carbon atoms, was thrust into the limelight of 
condensed-matter research since its experimental emer- 
gence in 2004^'^. Much of this attraction is due to its 
2D structure and contrary to any ordinary 2D material, 
having two Dirac cones in the Brillouin zone (BZ), where 
the conduction and valence bands touch. Charge car- 
riers with momenta near these two cones (known as K 
and K') have a unique linear energy dispersion and be- 
have like massless Dirac fermions. On the other hand, 
crystalline BLG^~^ has recently attracted a great deal of 
attention because of its unique tunable electronic prop- 
erties. It consists of two single-layer graphene (SLG) 
sheets separated by a small distance and can be produced 
by mechanical exfoliation of thin graphite or by thermal 
decomposition of silicon carbide. The low-energy quasi- 
particles in BLG behave as massive chiral fermions and 
are responsible for a plethora of interesting physics in- 
cluding broken-symmetry states at very weak magnetic 
fields when BLG is suspended to reduce disorder^ and 
anomalous exciton condensation in the quantum Hall 
regime^. Although the intrinsic BLG is a zero-gap semi- 
metal, it becomes a tunable band gap semiconductor^'^ 
when a gate voltage is applied. The band gap determines 
the threshold voltage and the on-off ratio of field-effect 
transistors and diodes, thereby making BLG more con- 
venient for applications in nano-electronic industry than 
SLG lO'ii. 

One of the fundamental problems of interest in 
graphene research is the indirect exchange interaction 
between two localized magnetic moments placed on this 



otherwise non-magnetic material. This carrier-mediated 
exchange interaction is known as RKKY interaction-'^^"-'^^ 
and it plays a significant role in the magnetic ordering 
of many electronic systems including spin glasses and al- 
loys. As it was originally studied for three-dimensional 
electron gas, it has also been studied for electron gas in 
one^^ and two^^ dimensions. Two main features of the 
long-range behavior of the interaction, measured by ex- 
change integral, J, for an electron gas is that it oscillates 
(in sign and magnitude) with the distance, between 
the moments, which exhibits ferromagnetic (FM) or anti- 
ferromagnetic (AFM) ordering and also decays^^'^^ with 
R. Both of these features have different functional forms 
depending on the dimension and generally, on the energy 
dispersion of the host material. For SLG, the RKKY 
interaction has extensively been studied For an un- 
doped SLG {Ep = 0) two main features are agreed upon, 
first, unlike an ordinary 2D metal with R~^ decay in the 
long-distance limit, J in undoped graphene falls off as 
R~^ and shows the 1 + cos[(X — K^) • i^]-type oscilla- 
tions with additional phase factors^^ depending on the 
direction of i^, and second, the moments on the same 
sublattice exhibit an FM interaction and an AFM cou- 
pling if placed on the opposite sublattices, as required 
by the particle-hole symmetry^^. The RKKY interaction 
for doped graphene shows a long-range behavior similar 
to that of ordinary 2D electron gas with another oscilla- 
tory factor emerging from the Dirac cones. It was shown 
that two characteristic momenta, kp and K — K' can be 
tuned to exhibit an unusual beating of the RKKY inter- 
action for certain magnetic moment arrangements^^. 

The RKKY interaction in BLG has also been addressed 
by several researchers^^'^^"^^. The local moment forma- 
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tion for adatoms on BLG using a mean-field theory of 
the Anderson impurity model has been studied by Killi 
et al^^. They showed that the RKKY interaction be- 
tween local moments can be varied by tuning the chemi- 
cal potential or by tuning the electric field as it induces 
changes in the band structure of BLG. The symmetry of 
the RKKY interaction on the bipartite lattice at half fill- 
ing has been discussed recently^^ and the distance depen- 
dence of the RKKY interaction has been briefly reported. 
Furthermore, Jiang et al^^, investigated the RKKY inter- 
action in multilayer graphene systems and they showed 
that the thickness of the multilayer influences the inter- 
action in a complicated manner and that the interaction 
couplings fall off as in long-range regime for BLG. 

However, the previous studies have only considered the 
RKKY interaction in the half-fllled {Ej, = 0) BLG. Con- 
sequently, signiflcant tunability feature of the RKKY in- 
teraction due to both doping, where the Fermi energy is 
no longer zero, and the perpendicular electric fleld, which 
gives rise to the gapped BLG, have not been addressed 
in the literature. Both of these cases are of paramount 
importance when it comes to possible spintronic applica- 
tions and will be the main focus in the present work. 

In this paper, we extend the Green's function (GF) 
technique used for the RKKY interaction in SLG^^'^^ to 
BLG. All cases of undoped, doped, unbiased and biased 
system are considered. We use the effective two-band 
Hamiltonian for the BLG^^^^ and for the flrst time, re- 
port the analytical expressions of the RKKY interaction 
for unbiased BLG in terms of the Meijer G-functions. 
We also present the numerical results of the interaction 
in the presence of a perpendicular electric fleld and show 
how the long-range behavior of the interaction can be 
tuned by the gate voltage. We explore a beating pattern 
of oscillations of the RKKY interaction in the four-band 
continuum model in which two conduction bands are par- 
tially occupied ( highly doped system). Furthermore, we 
discuss the discrepancy between the short-range RKKY 
interaction calculated from the two-band and that ob- 
tained from the four-band continuum model. 

The paper is organized as follows. In Sec. II, we in- 
troduce the general formalism and all the required GFs 
that will be used in calculating the RKKY interaction 
within the GF approach. In Sec. Ill and IV, we present 
the analytical results of the RKKY coupling within the 
two-band and four-band models, respectively. The main 
numerical results using two-band and four-band models 
and a brief account of the difference between the RKKY 
interaction within the four-band model and the effective 
two-band approximation are presented in Sec. V. Finally, 
we summarize the results in Sec. VI and draw some con- 
clusions. We have used the same GF method in studying 
the impurity states induced by a single vacancy in SLG, 
that includes the behavior of the a electrons as well as 
the TT electrons both using model and density-functional 
calculations. ^^'^^ 




FIG. 1. (Color online) Lattice structure of the BLG. The 
A (B) sublattices are indicated by red (green) circles with 
corresponding intra-layer and inter-layer hopping amplitudes. 
The bias voltage is denoted by V. 



II. THE MODEL HAMILTONIAN AND 
FORMALISM 

The BLG in Bernal stacking lattice shown in Fig. 1 
consists of two SLG lattices offset from each other in the 
xy plane with four atoms in the unit cell such that the top 
A-sublattice is directly above the bottom A-sublattice 
and it is between these pairs of atoms that the inter- 
layer dimer bonds are formed. The other two atoms do 
not have a counterpart on the other layer. We assume 
that the sp^-hybridized electrons of carbon atoms in each 
sheet are inert and only take into account the 2pz elec- 
trons which form the tt bands. 

We consider two magnetic impurities located at (a, 0) 
and R) and in contact interaction with the electrons 
of the biased BLG in Bernal stacking, where a and /3 
denote the sublattice indices (= Ai, 5i, ^2, ^2)- The 
tight-binding Hamiltonian of the system is given by 

= Ho + Hint 5 (1) 
where the Hamiltonian for the biased BLG, Ho, is given 

by 

1=1,2 i,a=A,B 

i j = l-3, 1 = 1,2 
-t^Yl^A^,R,^A2,R. - 73 X]4l,R,CB2,i^.+5,• 
-74 + c\,,r,cb2,r.+s,) + H{(2!] 

where / is the layer index, V is an external potential dif- 
ference between the layers, t = 2.9 eV is the intra-layer 
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nearest-neighbor hopping energy^^, the hopping energy 
between on-top sublattices in different layers is t_L = 0.3 
eV and furthermore, 73 = 0.12 eV denotes the hopping 
energy between not on-top sublattices between two lay- 
ers^. Another inter-layer second- nearest-neighbor hop- 
ping energies 74 = 0.04 eV and hence very small com- 
pared to t and can be ignored. The position vectors of 
three nearest neighbors of A-atom is denoted by Sj and 
a ^ 1.4 A is the carbon-carbon bond length. In the most 
general case, the on-site energies on the four atomic sites 
are no longer equal. They consist of independent pa- 
rameters to describe inter-layer asymmetry between the 
layers, an energy difference between two atoms in each 
layer and finally an energy difference between dimer and 
non-dimer sites. However, in this paper, we assume the 
equal on-site energies. The wave function can be written 
as a four-component spinor, V^Ai , V^Bi , V^As 7 V^Bs • this 
basis, the Hamiltonian of the biased BLG in Eq. (2) is 
represented as a 4 x 4 matrix given by 



Ho 



f V/2 f{k) 

f*{k) V/2 

t± Vif{k) 

\vif{k) v^rik) 



t± f4/*(fc)\ 

-v/2 r(fc) 

f(k) -v/2 J 



(3) 



where f{k) = -t J^Ui e 



ikSi 



V3 



73/t and V4 = 74/^- 



The interaction between the localized spins Si and ^2 
and the itinerant electron spins s is given by 



?^int = -A(Si- 81+52 -82). 



(4) 



Ignoring 73 and 74 in the Hamiltonian justifies the 
use of the BZ for SLG to describe the electrons mo- 
menta for the BLG. Therefore, we similarly describe the 
physics for those electrons with momenta in the prox- 
imity of the Dirac points Kjj = K^K'. In order to 
find the low-energy Hamiltonian near Dirac points, we 
expand the function f{k) using k = q -\- Kjj in pow- 
ers of q and keep only the linear term, which yields 
f{k) = f{q + Kd) — vpqe'^^^'^ where s = ± indicates 
the valley label, vp = 3ta/2 { h = 1 from now on) is the 

Fermi velocity of the electrons in SLG, q 



^Qx^Qy and 

Oq = tdiii~ ^ {qy / Qx) is the polar angle of q with respect to 
the X-axis chosen to be along the direction of X — K\ 
Furthermore, we consider vpq^V <C t± which allows us 
to eliminate the high-energy states perturbatively and 
simplify the Hamiltonian in Eq. (3) to two-band effec- 
tive Hamiltonian with states localized around Bi and B2 
sites. Doing so yields four bands, two degenerate bands 
in each of the two valleys K and K\ described by the 
effective two-band Hamiltonian^^ 



„ _ -1 / g2g-2is^,\ 











V/2 



,(5) 



where s = -hi for the K valley and s = — 1 for K' valley 
and m = t±/{2vp) and it is about 0.03me corresponding 
to a very small effective mass. The spinor is defined as 
= {ciBi^^^B2^ where creates an electron mostly at 



the Bi site with a small admixture from the other sites. 
We emphasize that the unperturbed Hamiltonians for the 
four-band and two-band models are denoted by T-Lq and 
iJo, respectively. 



A. The RKKY interaction Jcx,p{R) 

In the linear-response theory, the strength of the 
RKKY interaction, J, is found by two steps. First, using 
the Lippmann-Schwinger equation |^) = |^°)+G°V|^), 
one calculates the perturbed state |^), of the surround- 
ing electron gas (host material) at the unperturbed state 
1^^) due to the first moment. Si localized at the ori- 
gin and second, the first-order correction in the energy 
of this spin-polarized gas is found in the presence of the 
second moment, S2 localized at the lattice position 
viz., E{R) = (^|V(i^)|^). Here, = (£; + iO+ -Ho)"^ 
is the the unperturbed retarded GF. Therefore, the in- 
teraction energy may be written as 



E{R) = J{0,R)Si-S2, 



(6) 



with the RKKY interaction J(0, i^) being proportional 
to the static susceptibility, x(0, -R) viz.. 



j{o,R) = ^xio,R), 



(7) 



where the static susceptibility measures the proportional- 
ity between the perturbation 6V and the resulting change 
in the density Sn^ viz., xC^^'^^O = Sn{r) /SV{r^). 
It can be shown that written as^^ 

X(r,r') = -- Qm[G°iry,E)G°{r',r,Em 

J-00 

where G'>{r,r',E) = j:^Mr)r^{r'){E + iO+ - E^)-^ 
is the real-space matrix element of the retarded GF for 
a single spin channel with /i labelling the complete set 
of eigenstates of T-Lq. The factor 2 behind the integral 
counts for both spin channels. Eq. (8) is obtained 
by using the relationship between the charge density 
and the perturbed GF, viz., n(r) = J])!^^ IV^mI^)!^ = 
~f dE G{r^ r, E) and obtaining the charge dif- 
ference Sn{r) = n{r) — n^{r) induced by the perturba- 
tion 6V^{r^) from the approximated Dyson's equation 

The expression for the susceptibility in Eq. (8) 
can easily be extended to a system with several sub- 
lattice degrees of freedom, e.g. BLG. In a simi- 
lar definition for the magnetic susceptibility in the 
spin-density functional formalism, one can define the 
change in the density as^^ Sna^{r) = na^{r) — 
K/si'^) = T.a'(3' I (^^'Xa(3,a'(3'{r,r')Va'(3'{r') whcrc a 
or f3 denote the sublattice indices (e.g. ^41, 51, A2 
and B2 for BLG) satisfying the closure relationship 
/ |r, a)(r, a|(ir = 1 with the collective sublattice in- 
dex ly = ... and the perturbing potential is defined 
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as Vaf3{r^ r') = Vaf3{r)6{r — r'). Following similar steps, 
we find the generalized susceptibility as Xa/3,a'/3'{f^, '^^) = 
^J^:^dE ?sm[Gl^,iry,E)G%,ir',r,E)], and If we 
restrict the response only to the diagonal external poten- 
tial, the susceptibility in terms of the diagonal density 
matrix is given by Xq:/3(^5^') = ^^q:(^)/^^(^')5 which 
finally yields 



for sublattice components of the exchange integral as 

J«^(i^) = ^Xa/3(o,i^). (10) 

Knowing the real-space GFs, Eqs. (9) and (10) are the 
central formulas for calculating diff"erent sublattice com- 
ponents of the RKKY interaction in BLG. 



Xa0{r,r') = -- [ \E?sm[G%{r,r',E)G''^^{r',r,im 

^ J — OO 



Based on Eq. (9), for two magnetic moments one located 
at (q^, 0) and the other at (/3, R), we can re- write Eq. (7) 



B. Green's functions for the effective four-band 
Hamiltonian: Unbiased Case 

In the absence of the perpendicular electric field {V = 
0) , the unperturbed GF in momentum space correspond- 
ing to the four-band Hamiltonian T-Lq of Eq. (3) is repre- 
sented by 



1 

A 



-fik)nk)) 
^'-f{k)nk)) 

f{k)et^ 



e [e 



mr{k) 

-f{k)t^ 



f*{k)et^ \ 

nk)et^ r\k)t^ 

e - /(fc)/* (fe)) r (fe) - /(fc)/* (fc)) 

f{k) {e^ - f{k)r{k)) s {s^ - /(fe)r (fe) - ti)j 

(11) 



where e = E -\- iO~^ and A = 
(£2 - /(fc)r (fe) - etx) (s^ - /(fe)/*(fe) + 
zeros of which obtain the dispersion relation 
E{k) = ±t±/2 ± ^f{k)f'^{k) + (t±/2)2 leading to 
the celebrated Mexican- hat band structure of BLG. 
For sublattices, Ai, 5i, ^2, and the GF expres- 
sion in Eq. (11) has only six independent matrix 
elements, namely: G^^^^, G\^j^^, G%^b,^ G\^b^, 
^AiB2 ^BiSs' "^^^ corresponding matrix ele- 

ments of the real-space GF are obtained from the 
Fourier transformations of G^^(A;,£) elements, viz., 

G'i^{r,r',s) = Q^'^ Je'''<^--">G%ik,e)dk, where 
^BZ = (27r)^/l^ceii is the area of the first BZ with the 



area of the unit cell of SLG to be ftceW = 3A/3a^/2. 

Within the Dirac-cones approximation, the Fourier re- 
lationship simplifies and the general real-space GF ele- 
ment connecting the points (a, 0) and (P^R) is given by 

Gl^iO,R,e) =rtJdQ e-^''«[e-^«G0^(9 + K,£) 
+ e-'^'-^G%iq + K',e)], (12) 

from which the replacement R ^ —R yields the expres- 
sion for G^Q,(-R, 0, £:). The details of integrations for the 
Fourier integral in Eq. (12) are very similar to what has 
been reported in Ref. [24] and we will not repeat them 
here. Here, we just report the final results for the six 
matrix elements. They are 



G°AMO,R,e) = Ce(e-^^-« + e-^^'-^)[Ko{V^R) + Ko(V^R)] 

GXB,iO,R,e) = iC«F(e-^'^-^+'^« - e-''^'-''-"'^)[y^Ki{^/^R) + ^^K^{^^R)] 
GXA,(0,R,e) = Ce(e-^^-« + e-'^'-«)[Xo(\/^i?) -i^o(\/^i?)] 

GXb,{0,R,s) = -iCt;^ (e-^^-«-^^« - e-'^'-«+^^«)[v^ifi(\/^J?) - v^^i(v^^)] 
G%^sA^,R,e) = C(e-^^-« + e"^^' ■«)[(£ - t^)Koi^/^R) + (s + t^)Ko{V^R)] 

G%^B,{0,R,s) = ^ (e-^^--^-^^'^^ +e-^^'-«+2*''«)[a2i^2(V^i?) -/3'i^2(\/^i?)], (13) 



where Or = tan ^(y/x) is the polar angle of the direction direction of K — K\ ( = —TrVp^fl^^, o? = v^[e^ —et\})^ 
of R with respect to the x-axis chosen to be along the 
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0- 



2 _ 



^(^^ + £t^) and K^{x) is the modified Bessel where the function C^jp in the bracket is the Meijer G 



function of the second kind and order of // = 0, 1, 2. 



function^^ . 



C. Green's functions for the effective two-band 
Hamiltonian 

Using the Hamiltonian in Eq. (5), the momentum- 
space matrix representation of the retarded GF defined 
as G^{e) = {s — Ho)~^, is given by 



s + V/2 



-V/2 



(14) 



£^ with the band energy dispersion 



where A' = a 

Sq = zb^ + It should be noticed that the re- 
tarded GF in Eq. (14) has two independent terms, 
namely G%^j^^ {q, e) and G%^j^^ {q, e) since G%^j^^ {q, e) = 
G%^^^{q^e). Furthermore, G%^^^{q,s) can be obtained 
from G^_^Q_^{q^e) by replacing V —V. 

Points on the same layer- Similar to the previous sec- 
tion, we use Eq. (12) to find the corresponding real- 
space GFs. The Fourier integral can be evaluated in two 
ways. As the first method, we can plug the expression 
for G%^Q^{q,s) given in Eq. (14) into Eq. (12) and in- 
tegrate. After some algebra, we finally obtain the GF in 
terms of the Meijer G- function as 



27r 



2VV 



(2e+V) ^3 



5,0/^ 1 1 ^ 
4 ^^U, 2, 25 ^ 



iK R ^ 

rn^R^/T/2 

256 



-4.2))] ,(15) 

I 



As for the second method, we re- write the cor- 
responding momentum- space GF as b i^^^) ~ 
m2(2e + F)Ei=i,2rn^ + (-l)V)-' where ^ = 
V^rn'^s'^ — m^V^. After the integrations, the GF reads 



'^^^f^ [K,{V~m - KoiV^R)] , (16) 



We emphasize that ^ is a c-number depending on the 
range of the energy. Although both of the expressions in 
Eqs. (15) and (16) are equivalent, we will be using the 
one in terms of the modified Bessel function, which will 
be practical for numerical analysis. 



Points on different layers- 



In this case, we write 
__.2V-iy(e + (-l)V)-^. Fol- 
lowing the steps elaborated for the same layer, we obtain 
two equivalent expressions for the real-space GF corre- 
sponding to this case in terms of both Meijer G- and 
modified Bessel functions as 



-i{K'-R-2eR) 



3,0/^ 1 - 
0,4 \ ^-> 2 ' - 



256 



2m7r 



i{K-R+2eR) 



i{K'-R-2eR) 



[K2[VIR) + K2{^R)\ 



(17) 



III. RKKY INTERACTION FROM THE 
TWO-BAND MODEL 

In this section, for the moments on the same 
layer, we use the expressions for the real-space GFs, 
G^,B,(0,il,e,F) and G^^B^(il,0,£,y) from Eq. (16) 
and then obtain the RKKY interaction JbibAR) using 
Eqs. (9) and (10). 



A. Moments on the same layer: JbiBi{R) 

General case- Using Eq. (16) and substituting both 
G%^^^{0,R,e,V) and G^^^^ (i^, 0, 5, V) in Eq. (9), the 
corresponding susceptibility reads 

— 167r?T7^ 

XB,B,(0,il) = — ^2 ^b^bMV,R,Ef) (18) 



where ^BiBi = 1 + cos [{K — K') - R] and the integral 
Ii is given by 

[KoiV^R) - KoiV^R)]' ■ (19) 

Ii may not be analytically evaluated for arbitrary val- 
ues of the gate voltage, V and Fermi energy Ep- How- 
ever, for the special case of unbiased system, ^ = one 
can manage to find the analytical expression for Ii as 
elaborated in the following subsection. 

Special case: Unbiased BLG- Now, we split the inte- 
gral in Eq. (19) into two parts, viz., /^^ = /^^ + j^"" , 
where the first term accounts for the valance electrons 
(undoped case) and the second for the conduction elec- 
trons (doped case). Let's denote the first integral by 
/q. Special care must be taken to consider the prin- 
cipal value of the complex square roots of the com- 
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plex variables. For instance, ^ = 2m\J = ±2m6 for 
^ > and ^ < 0, respectively. Then, we introduce 
new variables y = -^^mEB? accordingly for both inte- 
grals such that 2/ > in each and express the modified 
Bessel function with complex argument in terms of the 
Hankel functions and then Bessel functions of first and 
second kind using K^^z) = I^U^y^ {ze'"" 1^) = 

-2-ii7re-^^^/2iji^^(ze-^^/2) with hI'^'^\z) = J^{z) ± 
i Yiy{z). In particular, we use Ko{y/y =t iO+) = Ko{-y/y) 
and Ko{V-y±m) = {-7r/2) {Yo{^) ± iM^)) . 
After some algebra, /q simplifies to 

I ir dyMVy)yo{Vy)]- (20) 

Both integrals in Eq. (20) are diverging; however, af- 



ter using some regulatory cut-off functions^^ they are 
evaluated to one and zero, respectively, which yields 
Jo = Ti{2mR^)~^ . Plugging this result into Eqs. (18) 
and (7) immediately gives the RKKY interaction for the 
unbiased and undoped BLG as 



1 -h COS [{K - K') ' R] 



(21) 



where C = 3A^(167r^t^) ^t_L is a positive parameter, 
which means that Jb^bS-^) represents an FM interac- 
tion between the moments. The power-law R~'^ decay of 
the RKKY interaction in Eq. (21) clearly shows that the 
undoped and unbiased BLG behaves like an ordinary 2D 
electron gas, the result that have also been reported in 
other studies^^'^^-^o. 

As for the general doped case, after similar simplifica- 
tions, we obtain 



h{V = Q,R,XF) = 



2mi?2 



Jo 



dyM^)KQ{^) 



2 Jo 



dyM^)Yo{y/y) 



(22) 



where xp = 2mEFR^ = kj^R'^. Both inte- 
grals in Eq. (22) can be expressed^^ in terms of 
the Meijer G-functions, viz., J^^ dyJo{-^)Ko{y/y) = 



0,0 



1 

' 2' 2' 







xl/64 



and similarly. 



Xf 



We find the asymptotic expansion of the Meijer G- 
functions^^ in Eq. (22) and eventually obtain the RKKY 
interaction for the large distances as 



lim UV = ^,R,kF) = ^^-^ 



V^e-^-^ cos (kFR) + I sin (2A:^it!) - J^^^^) 
2 SkFR 



(23) 



The peculiar feature of the interaction at this limit is 
the exponential decay along with the power-law decay. 
However, the exponential term does not survive for very 
large distances. The reason for such an exponential decay 
can probably be explained based on the chiral character- 
istics of the quasiparticle in BLG and the fact that the 
forward scattering is forbidden in the system. 



B. Moments on different layers: Jb-iB^^R) 

General case- In this section, we consider the situation 
in which the magnetic moments are located on different 
layers. Using Eqs. (17) and (9), the corresponding sus- 
ceptibility, XB^B^i^^R) is given by 



— 167rm^ 

XBiB.(0,i?) = $B,B, I2{V,R,Ef), (24) 



where ^BiBs = 1 + cos [{K - K') • R + A6ii\ and the re- 
maining integral is given by 



h{V,R,EF) 



i-Ef 
I ' 

J —00 



dE 



K2{^/^R) ^ K2{^/^R)h5) 



Similar to the case of the moments on the same layer, we 
can find the analytical expression for the RKKY interac- 
tion for the unbiased BLG. 

Special case: Unbiased BLG- Here, again we split 
the integral of Eq. (25) into the undoped and 
doped parts and perform the same type of calcu- 
lation s as ex plained before. In parti cular, we use 
i^2(vV±^) = K2{^) and K2{V-y±iO+) = 
(7r/2) (Y2{-s/y) ^ iJ2{-\/y)) ' Therefore, the integral cor- 
responding to the undoped part is denoted by /q and it 
reads 



2mi?2 



[!^dyJ2{^)K2iVy) + 



(26) 
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Resorting to cut-oflF function scheme, both diverging inte- 
grals in Eq. (26) give one and — 4/7r, respectively, which 
yields = -7r(2mi?2)-i Finally, using /q in Eqs. (24) 
and (7) we obtain the RKKY interaction for the unbiased 
and undoped BLG for the moments of different layers as 

Jb.b.{R) = C . (27) 

As C > 0, Jb^bS-^) signifies an AFM interaction be- 
tween the moments on different layers. 

The comparison between Jb^B2^-^) Jb^bS-^) 
given in Eq. (21) reveals a very subtle point. Apart 
from their different oscillatory Dirac-cones factor, which 
are both bound and positive, Jbibi and JbiB2 ^tie 
same magnitude C and opposite sign. We recall that 
for the case of undoped SLG, the RKKY interaction for 
the moments on the opposite sublattices is AFM and its 



magnitude is three times larger than that of for the same 
sublattice, namely. Jab = ^Jaa^^'^^'^^: which means 
the AFM ordering is more favored for SLG. It was rea- 
soned by Saremi^^ that this commensurate feature of the 
RKKY interaction must be the case for any bipartite sys- 
tem with particle-hole symmetry. By analogy, our results 
for JbiBi and JbiB2 be interpreted as the signature 
of the bipartite nature of the system and the particle- hole 
symmetry present in the effective two-band Hamiltonian 
Hq in Eq. (5). Although it appears that the unbiased 
BLG bears the same symmetry, an attempt to prove the 
theorem particularly for this system will be insightful. 
To our knowledge such proof does not exist. 

Similarly, we can find the analytical expressions of the 
interaction for the doped case. Following same steps as 
discussed previously, Eq. (25) simplifies to 



2mE? 



j^J dyj2{Vy)K2{Vy) + f dyj2{Vy)Y2{Vy) 



(28) 



where the first and second integrals in Eq. (28) 
are evaluated as j^^ dyJ2{^)K2{y/y) 



2,1 



1 5 

1,3,-1,0 



Xf 



-1/2 



3,1 



xpGt: 



5 I 1 1 



x|/64 



and 



respectively. The long-distance expression is given by 



TT 



lim hiV = Q,R,kF) = -. „2 



V2e-'^^ cos {UfR) - l sin {2fc^i?) - ^'^/'t^"^^ " 
2 okpR 



(29) 



Comparing Eq. (23) and (29), we note that the long- 
range functional form of the RKKY interaction in unbi- 
ased BLG for the moments on different layers is the same 
as for those on the same layer. 



four-band model. For the biased case the same analysis 
can be made, which will not be presented here. 



IV. RKKY INTERACTION FROM FOUR-BAND 
CONTINUUM MODEL: UNBIASED BLG 

In this section, we report the expressions for the RKKY Plugging the GFs from Eq. (13) into Eq. (9), the 
interaction for the unbiased and doped BLG using the corresponding susceptibilities are given by 



J — c 

XAiBi(2)(0:^) =^^AiBi(2) / 



2c^ 

X \sm 



Koi\l-x^+*-^x±Ko{\l-x 



t±R 
vf 



x) 



dx 



t±R 



x'^ H xKi{\ —x'^ H x) dz \ —X 



Vf 



t±R 



xK,{J~x^ 



t_i_R 

VF 



X) 



dx 
(30) 
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XBiBi(0,-R) = A$BiBi [ 

J — ( 



dx (31) 



J — oo 



2 tj^-R 



t'F V '^F 



dx, 
(32) 



where x RE/vy, A = -AnVt^^v^^ 

and = REy/vy, ^a^a^ = ^BiBi = 

= 1 + cos[(K-ii:o = 

1 - cos [(K -K') 20 r] and ^^iS^ 

l + cos[(X-XO •i^ + 4(9i^]. 

To calculate the RKKY interaction in the four-band 
continuum model, Eqs. (30)- (32) may be evaluated nu- 
merically. The results of Eqs. (31)- (32) and those ob- 
tained from the two-band model are compared in the 
following section. 

V. NUMERICAL RESULT 

In this section, we present our main calculations for 
the exchange coupling of the RKKY interaction evaluat- 
ing Eqs. (19), (25) and (30)-(32). The general features of 
the exchange coupling, basically the dependence of the 
RKKY interaction on the distance R have been numeri- 
cally studied previously^^ for the unbiased and undoped 
BLG. We present, on the other hand, our numerical re- 
sults of Ii{V, R, Ey) and hiV, R, Ey) for biased BLG in 
two different interesting regimes namely doped and un- 
doped graphene where Ey = ^ and 7^ 0, respectively. 
We provide a comparison between the results obtained 
within the four-band and the two-band continuum mod- 
els in unbiased BLG systems and discuss the discrepancy 
between two models. 



A. Unbiased and doped BLG 

In previous Sections, we found the analytical expres- 
sions of the RKKY interaction for the unbiased BLG 
within the two-band model. We showed that regardless 
of the Dirac-cones oscillatory term represented by ^a,(3i 
the main difference between the interactions for the mo- 
ments on the same and different layers is due to their 
sign, which results in FM interaction between impurities 
on the same layer BiBi and AFM interaction in B1B2 
case. By doping BLG, the strength of the RKKY inter- 
action decreases and a new oscillatory behavior starts. 
Therefore, the RKKY interaction changes sign as a func- 
tion of distance. 

Fig. 2(a) shows IiiV = O^R^Ey) integral as a func- 
tion of the distance, R for different doping values. The 
period of oscillation and the speed of a decay depends 



strongly on the Fermi Energy. For non-zero Ey, the inte- 
gral Ii shows a quite different behavior as it exhibits an 
oscillatory behavior as a function of R with decreasing 
amplitude and a period given by tt/Zcf- We compare the 
analytical results of Eq. (23), plotted as solid lines, with 
the numerical evaluation of Eq. (22), plotted as symbols, 
to show their difference at short distance while reaching 
each other quite well at large R regions. A comparison 
between i?-dependence of the integral Ii and that of I2 
for Ey = 0.05 eV in Fig. 2(b), shows their difference at 
short distance while reaching each other approximately 
as R increases. Similar to /i, at finite Ey, the integral 
I2 has an oscillatory behavior as a function of R, with 
decreasing amplitude and a period given by tt/Zcf- 

Fig. 3 shows a comparison between the RKKY inter- 
action that calculated by the four-band model given by 
Eqs. (31), (32) and those obtained by two-band model. 
For an undoped case, it can be seen from Fig. 3(a) 
that two results are matched for FM interaction, (/i) 
while there is a discrepancy between results at short- 
range regime which is controlled by value of t± for AFM 
interaction, {I2). This result is very pertinent to the 
conclusion stated in Ref. [40] where the authors show 
that based on the charge-charge response function calcu- 
lations, the density-sum and density- difference fluctua- 
tions in BLG crossover from those of an unusual massive- 
chiral single-layer system to those of a weakly coupled 
bilayer as carrier density, wave vector, and energy in- 
crease. Fig. 3(b) shows the Ii for unbiased and doped 
BLG when Ey = 0.05 eV calculated by the four and two 
band continuum models. It is clear that the period of the 
oscillation in the four-band model is different from that 
of the two-band model. This is because for a certain and 
small Fermi energy value, the associated Fermi momen- 
tum is larger than the value obtained in the two-band 
model. Note that the electronic dispersion relations in 
the four-band model (roots of A defined in Eq. (11)) 
may be written as £^i,2(^) = \/'^y^'^ ^-t\/4:± and 
E^^^{k) = y^v^'^ ^t\/4±t^/2. Therefore, the period 
of the oscillations depends on the model. 

In the four-band model, depending on the doping level, 
the Fermi energy can have either one or two intersections 
with the conduction-band Fermi surfaces. By increasing 
the Fermi energy, the proper results of the quasiparticle 
excitation are captured by the four-band model. This 
point is demonstrated in Fig. 3(c) where we show the re- 
sults of the I2 for Ey = 1 eV, for which the Fermi energy 
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FIG. 2. (Color online) (a) Ii{V = 0,R,Ef) as a function of 
the distance R for different doping values obtained for the 
two-band model and unbiased BLG using Eq. (22). The 
function R^Ii{V — 0,R,Ef = 0) is a constant in agreement 
with that result obtained in Ref . [26] . Solid lines refer to the 
analytical results of the asymptotic behavior from Eq. (23) 
which are compared to the numerical evaluation of Eq. (22), 
plotted as symbols, show their difference at short distances 
while reaching each other as R increases, (b) The strength 
of the interaction for both BiBi [Eq. (22)] and B1B2 [Eq. 
(28)] for undoped and doped systems for = 0.05 eV as a 
function of the distance. 



intersects two conduction bands. It is clear that in this 
case the results obtained by full band are completely dif- 
ferent from those calculated by the low-energy excitation 
method. 

The intersection of the Fermi energy with the bands 
denoted by £^1(3) creates two Fermi surfaces. Because 
the RKKY interaction is fundamentally determined by 
the geometrical features of the Fermi surface of the 
host material, a somewhat more complicated behavior 
of the RKKY coupling for a highly doped BLG can oc- 
cur. As the result, we observe that oscillations of XBiBi 
exhibit a beating pattern with two characteristic peri- 
ods associated with the two Fermi momenta defined as 
/cfi(2) = V-^F^ EFt±/vF' Fig. 4 shows this beating pat- 




FIG. 3. (Color online) (a) Comparison between the results 
of the RKKY interaction that obtained by the two-band 
model [Eq. (28)] and that calculated by the four-band model 
[Eq. (32)] for I2 at V = Ef = 0. There is a discrepancy 
between two approaches which basically comes from the off- 
diagonal inter-layer tunneling term t±. In the inset, a com- 
parison of results between two theories for the case of Ii . (b) 
/i(y = 0,R,Ef) for = 0.05 eV and (c) hiV = 0,R,Ef) 
for Ef = 1 eV as a function of the distance R obtained by the 
two-band model [Eqs. (22) and (28)] and that calculated by 
the four-band model [Eqs. (31) and (32)], respectively. The 
proper results of the quasiparticle excitation are captured by 
the four-band model, by increasing the Fermi energy. 
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360 



FIG. 4. (Color online) The susceptibility XBiBi as a func- 
tion of the distance between impurities on the same sublattice 
along armchair direction obtained from the four-band model 
[Eqs. (31)]. The existence of two different periods (beating 
pattern) in doped BLG for certain values of = 1 eV is 
clear in this figure. 
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FIG. 5. (Color online) (a) The integral Ii{V,R,Ef = 0) as 
a function of R for different gate voltages. The function falls 
off rapidly and oscillates slightly for finite V values, (b) the 
same as (a) but for the integral hiV, Ef = 0). 
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FIG. 6. (Color online) (a) Ii{V, R, Ef) as a function of R 
for different gate voltage at Ef = 0.1 eV. To emphasis the 
amplitude value of Ii(V, R, Ef), R^ Ii(V, R, Ef) is shown for 
different gate voltages, (b) The same as (a) for hiV, R, Ef). 
In the inset of (b), solid lines refer to the analytical results 
of Eq. (33) and are compared to the numerical evaluation of 
Eq. (25), plotted as symbols, show their difference at short 
distance while reaching each other quite well as R increases. 
Here, a = 0.42, /3 = 1.6 and a = 0.33, /3 = 1.45 for V = 0.1 
and 0.14 eV, respectively. 



tern the RKKY interaction as a function of the distance 
between impurities. 



B. Biased BLG, doped and undoped 

By turning on the gate voltage perpendicular to the 
system, a new type of dispersion relation of the band is 
emerged. At zero Fermi energy, due to a gap opening and 
consequently removing the available energy states for the 
mediating electrons, the response function between elec- 
trons decreases, and as a result, the RKKY interaction 
decreases much faster than R~^. 

Fig. 5(a) shows the integral /i(F, i?, = 0) as a func- 
tion of R [Eq. (19)] for different gate voltages in the 
two-band continuum model. The function Ii{V^R^Ef = 
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0) decreases by increasing the bias voltage and oscil- 
lates slightly around its zero value. Decaying structures 
and oscillations depend on the bias voltage. Similarly, 
Fig. 5(b) shows the integral /2(V,i^, E'f = 0) as a func- 
tion of R [Eq. (25)] for different gate voltages. The 
function decays as R increases and remains negative. 

One goal of the present work is to understand the 
RKKY interaction in a doped BLG system. For this pur- 
pose, we consider a finite gate voltage together with the 
finite Fermi energy to calculate the RKKY interaction. 
Fig. 6 shows the integrals Ii{V,R,Ef) and hiV.R.Ep) 
as a function of R for different gate voltages at given 
= 0.1 eV. Our results show that Ii{V^ R^ E^) is sen- 
sitive to V and by growing it around 2Ep values, the 
amplitude and also the wavelength of the oscillation of 
/i(y,i?, £^f) increases. However, /2(V,i?, ^f) slightly 
changes with the gate voltage. Similar to the case of 
unbiased and doped BLG, the integral Ii for non-zero V 
exhibits an oscillatory behavior as a function of R with 
a period now controlled by both the Fermi momentum 



■a 2mR'^ 



where ky = [{2m E^)^ — m'^V^Y ^, /3 and a are param- 
eters controlled by ^f and V. In the inset of Fig. 6(b), 
solid lines refer to the analytical results of Eq. (33) com- 
pared to the numerical evaluation of Eq. (25), plotted as 
symbols, showing their difference at short distance while 
reaching each other quite well as R increases. 



VI. SUMMARY 

We have studied the effect of the bias voltage on the 
RKKY interaction in doped and gapped BLG. Our ap- 
proach is based on the lattice Green's function tech- 
nique. Near the Dirac points, charge carriers in BLG 
have parabolic energy spectrum with a finite density of 
states at zero energy, similar to the conventional non- 
relativistic electrons. On the other hand, these quasipar- 
ticles are also chiral and described by spinor wave func- 
tions. Therefore, the dependence of the RKKY interac- 
tion on the position vector R between two local magnetic 
moments is not only controlled by the dispersion rela- 
tion but also by the chirality, which makes it directional- 
dependent as also shown to be the case for SLG ^^'^^ by 
the phase factors ^a^- 

Similar to SLG, we report the ferromagnetic in- 
teraction for moments on the same layers and anti- 
ferromagnetic coupling for those placed on the opposite 
layers in unbiased and undoped BLG. We associate this 
feature to the particle-hole symmetry and the bipartite 
nature of the lattice within the two-band model as argued 



and the gate voltage. Our numerical results show that 
the period of the oscillations can be fitted quite well by 
ir/ky where ky = ((2m£^F)^ — ^^^^)''"^^- One interesting 
feature in this case is that the long-range behavior of the 
RKKY interaction for the impurities on the same layer 
is similar that of a standard 2D electron gas. Another 
interesting feature is the enhancement of the RKKY in- 
teraction by increasing the gate voltage illustrated by 
{R^h) results in the inset of Fig. 6(a). Since the RKKY 
interaction decays rapidly, it is almost difficult to mea- 
sure it experimentally. Based on our results, here we 
proposes that the tuning of the gate voltage to a certain 
value, 2Ey-, will noticeably enhances the strength of the 
RKKY interaction and thereby makes it accessible for 
experimental probes. 

Finally, we find that our numerical results for large dis- 
tances between two impurities located on different layers 
[Eq. (25)] can be faithfully fitted by an analytical expres- 
sion very similar to the unbiased case given in Eq. (29). 
This asymptotic fit is given by 



(33) 

I 

in Ref. [19]. 

For the unbiased and doped case, we managed to find 
the analytical expressions of the RKKY interaction in 
terms of the Meijer G-functions and their long-range 
behavior was also reported. The salient feature of the 
asymptotic behavior is that the power-law decay R~^ 
is accompanied by an exponential factor as JbiBi^2) ^ 
TR~^ cos{kFR)[e-^^^ ± 2-1/2 sin(A:Fi^)]. It was shown 
that the mediating carriers of a gapped graphene^^'^^ or 
SLG with disorder'^i produce an exponential decay in 
the RKKY interaction; however, for a pristine unbiased 
BLG, which is gapless, we associate this exponential de- 
cay to the chiral nature of the carriers in the system. 

We have supplemented the results from the two-band 
model for the unbiased case with our calculations using 
the four-band model to identify the validity of the two- 
band model and the discrepancy between both models. 
Within the four-band model, when the system is highly 
doped, the application of the two-band model is ques- 
tionable. In this regime, the main features of the RKKY 
interaction are only captured by the four-band model. In 
low-energy region, we have shown that the two models are 
different at the short-range of the distance between im- 
purities located on different layers. In addition, we have 
observed that the oscillations of XBiBi exhibit a beating 
pattern with two characteristic periods associated to the 
two Fermi momenta. 

For the biased and doped BLG, we have shown that 
the gate voltage and the Fermi energy can vary indepen- 
dently to determine the RKKY interaction. One of the 



V2e 



-kvR 



COS (kyR) — asin{2kYR) — (3 



cos {2kvRy 
k^ 



12 



References 


System 


RKKY interaction tor same sublattice 


RKKY interaction tor ditterent subtattices 


Ref. [16] 


2DEG 


R sm{2kFR) 




Ref. [43] 


2DEG+ impurity 


R sm.{2kYR)e 




Ref. [24] 


SLG (^F = 0) 


-R~^^AA 


3i?~^$AB 


Ref. [25] 


slg(£;f / 0) 


-R-'^ sm{2kF R)^AA 


R-'^ sm{2kF R)^AB 


Ref. [26,30] 


BLG(^F = 0,V = 0) 


-R-'^AA 


R-'^AB 


Present work 


BLG(^F 7^ 0, y = 0) 


-R-^ cos{kFR)[e-''^^ + 2-^/2 sm{kFR)]^B,B^ 


R-^ cos(kFR)[e-''^^ - 2"^/^ sm{kFR)]^B^B2 


Present work 


blg(£;f / 0, y / 0) 


-R-^sm{2l)^B^B^ 


R-^[e-^cos{l) asm{2l) ]$5iB2 



TABLE 1. A breakdown of ttie results on ttie scaling form of the RKKY interactions in two dimensional electron gas ( 2DEG), 
SLG and BLG. The RKKY interactions are proportional to values given in the third and fourth columns, a and /3 are parameters 
controlled by ^f and V. The parameter / = kyR where ky = ((2m^F)^ — m^V^)^^^. The functions ^aa and ^ab are given 
by 1 + cos[(ii: - K') • R] and 1 + cos[(ii: - K') • i^ + tt - 20r], respectively. 



fascinating features in this case is the possibility for the 
enhancement of the interaction by tuning the gate volt- 
age and/or the Fermi energy, which opens an avenue to 
probe the interaction experimentally. We have obtained 
the asymptotic behavior of the RKKY interaction ana- 
lytically for each case and the expressions are given in 
Table 1. In order to compare the RKKY interaction in 
BLG with an ordinary 2DEG we have reported the inter- 
action in clean 2DEG^^ and in the presence of disorder 
where the exponential decays is introduced^^. 



VII. ACKNOWLEDGMENTS 



MS thanks Nico Temme for useful discussions on the 
Meijer G-functions. The work at the University of Mis- 
souri was supported by the U. S. Department of Energy 
through Grant No. DE-FG02-00ER45818. 



asgari@ipm.ir 

^ K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. 

Zhang, S. V. Dubonos, L V. Grigorieva, and A. A. Firsov, 

Science, 306, 666 (2004). 
2 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. 

Novoselov and A. K. Geim, Rev. Mod. Phys. 81, 109 

(2009). 

^ C. Berger, Z. Song, T. Li, X. Li, A. Y. Ogbazghi, R. Feng, 

Z. Dai, A. N. Marchenkov, E. H. Corad and P. N. First, J. 

Phys. Chem. B 108, 19912 (2004) . 
^ K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. 

L Katsnelson, L V. Grigorieva, S. V. Duboson and A. A. 

Firsov, Nature (London) 438, 197 (2005) . 
^ Y. Zhang, Y.-W. Tan, H. L. Stormer and P. Kim, Nature 

(London) 438, 201 (2005) . 
^ O. Vafek, Phs. Rev. B 82, 205106 (2010); F. Zhang, H. 

Min, M. Pohni and A. H. MacDonald, Phys. Rev. B 81, 

041402 (2010) . 
^ K. S. Novoselov, E. McCann, S. V. Morozov, V. L Fafko, 

M. L Katsnelson, U. Zeitlet, D. Jiang, F. Schedin and A. 

K. Geim, Nature physics 2, 177 (2006); E. McCann and V. 

L Falko, Phys. Rev. Lett. 96, 086805 (2006); Ya-Fen Hsu, 

Guang-Yu Guo, Phys. Rev. Lett. 82, 165404 (2010) . 
^ Kin Fai Mak, Chun Hung Lui, Jie Shan, Tony F. Heinz, 

Phys. Rev. Lett.102, 256405(2009) . 
^ A. B. Kuzmenko, et al. Phys. Rev. B 80, 165406 (2009) . 

Taisuke Ohta, Aaron Bostwick, Thomas Seyller, Karsten 

Horn, Eh Rotenberg, Science 313, 951 (2006); Eduardo V. 

Castro, K. S. Novoselov, S. V. Morozov, N. M. R. Peres, 

J. M. B. Lopes dos Santos, Johan Nilsson, F. Guinea, A. 

K. Geim, A. H. Castro Neto, Phys. Rev. Lett. 99, 216802 



(2007) ; J.B. Gostinga, H.B. Heersche, X. Liu, A.F. Mor- 
purgo, L.M.K. Vandersypen, Nature Mater. 7, 151 (2008); 
W. Yao, D. Xiao, and Q. Niu, Phys. Rev. B 77, 235406 

(2008) . 

B. R. K. Nanda and S. Satpathy, Phys. Rev. B 80, 165 430 

(2009) . 

M. A. Ruderman and C. Kittel, Phys. Rev. 96, 99 (1954). 
T. Kasuya, Prog. Theor. Phys. 16, 45 (1956). 
K. Yosida, Phys. Rev. 106, 893 (1957). 
Y. Yafet, Phys. Rev. B 36, 3948 (1987) 
B. Fischer and M. W. Klein, Phys. Rev. B 11, 2025 (1975). 
^'^ M. A. H. Vozmediano, M. P. Lopez-Sancho, T. Stauber 
and F. Guinea, Phys. Rev. B 72, 155 121 (2005). 
V. K. Dugaev, V. L Litvinov and J. Barnas Phys. Rev. B 
74, 224 438 (2006). 

S. Saremi, Phys. Rev. B 76, 184430 (2007). 

L. Brey, H. A. Fertig and S. Das Sarma, Phys. Rev. Lett. 

99, 116802 (2007). 

J. E. Bunder and H.-H. Lin, Phys. Rev. B 80, 153414 
(2009). 

2^ A.M. Black-Schaffer, Phys. Rev. B 81, 205416 (2010). 
F. Parhizgar, R. Asgari, S. H. Abedinpour and M. Zareyan, 
arXiv:1211.2013 

M. Sherafati and S. Satpathy, Phys. Rev. B 83, 165425 
(2011). 

M. Sherafati and S. Satpathy, Phys. Rev. B 84, 125416 
(2011). 

2^ E. Kogan, Phys. Rev. B 84, 115119 (2011). 

M. Sherafati and S. Satpathy, AIP Conference Proceedings 
1461, 24 (2012). 



13 



E. H. Hwang and S. Das Sarma, Phys. Rev. Lett. 101, 

156802 (2008). 

M. Killi, D. Heidarian, and A. Paramekanti, New J. Phy. 
13, 053043 (2011). 

L. Jiang, X. Li, W. Gao, G. Yu, Z. Liu, and Y. Zheng, J. 
Phy. C 24, 206003 (2012). 

B. R. K. Nanda, M. Sherafati, Z. S. Popovic and S. Satpa- 
thy, New J. Phys. 14, 083004 (2012). 

M. Sherafati and S. Satpathy, Phys. Status Sohdi B 248, 
2056 (2011). 

E. McCann and V. I. Falko, Phys. Rev. Lett. 96, 086805 
(2006). 

J. Nilsson, A. H. Castro Neto, F. Guinea, and N. M. R. 
Peres, Phys. Rev. B 78, 045405 (2008). 
Eduardo V. Castro, K. S. Novoselov, S. V. Morozov, N. 
M. R. Peres, J. M. B. Lopes dos Santos, Johan Nilsson, F. 
Guinea, A. K. Geim, and A. H. Castro Neto, J. Pjys. C 
22, 175503 (2010) . 



L. M. Malard, J. Nilsson, D. C. Elias, J. C. Brant, F. 
Plentz, E. S. Alves, A. H. Castro Neto, and M. A. Pimenta, 
Phys. Rev. B 76, 201401 (2007). 

P. Mohn, Magnetism in the Solid State: An Introduction 
(Springer, 2005) . 

J. Kiibler, Theory of Itinerant Electron Magnetism (Ox- 
ford University Press, USA, 2009) . 

C. S. Meijer, Nederl. Akad. Wetensch. Proc. Ser. A 49, 
344 (1946); Y. L. Luke, The special functions and their 
approximations (Academic Press, 1969), vol. 1, p. 143, 171, 
191 and 233; J. Fields, Math. Comput. 26, 757 (1972). 

G. Borghi, M. Pohni, R. Asgari and A.H. MacDonald, 
Phys. Rev. B (R) 80, (2009) . 

H. Lee, J. Kim, E.R. Mucciolo, G. Bouzerar and S. Kette- 
mann Phys. Rev. B 85, 0754250 (2012). 

E. Kogan, Graphene 2, 8 (2013). 

U. Larsen , J. Phys. F. 15,101 (1985) . 



